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ABSTRACT 

Using a three-dimensional magnetohydrodynamic model, we simulate the 
magnetic reconnection in a single current sheet. We assume a hnite guide held, 
a random perturbation on the velocity held and uniform resistivity. 

Our model enhances the reconnection rate relative to the classical Sweet- 
Parker model in the same conhguration. The efficiency of magnetic energy con¬ 
version is increased by interactions between the multiple tearing layers coexisting 
in the global current sheet. This interaction, which forms a positive-feedback 
system, arises from coupling of the inhow and outhow regions in diherent layers 
across the current sheet. The coupling accelerates the elementary reconnection 
events, thereby enhancing the global reconnection rate. The reconnection estab¬ 
lishes hux tubes along each tearing layer. Slow-mode shocks gradually form along 
the outer boundaries of these tubes, further accelerating the magnetic energy con¬ 
version. Such positive-feedback system is absent in two-dimensional simulation, 
three-dimensional reconnection without a guide held and a reconnection under a 
single perturbation mode. We refer to our model as the “shock-evoking positive- 
feedback” model. 

Subject headings: (magnetohydrodynamics;) MHD, plasmas, magnetic reconnection, - 
methods:numerical 
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Introduction 


Magnetic reconnection is considered to sonrce the rapid conversion of magnetic energy 
in varions solar coronal activities with extremely high Lnndquist nnmber S = Lva/^I ~ 10^^, 
where L, va and fj denote the cnrrent sheet length, Alfven speed and magnetic di ffnsivit y, 
respectively. The two classical reconnectio n models are t he Sweet-Parker model f Sweetl 


19581 : 


Parker 


19631 ) and Petschek’s model flPetschekl 119641) . In the Sweet-Parker model, 


the reconnection rate scales as Vinfiow/^A ~ several orders of magnitnde slower 

than reqnired in many solar and astronomical applications. The Petschek’s model yields a 
snfficiently fast reconnection rate (~ I/Ins') by virtue of the localized diffusion region and 
the extended slow-mode shocks. The localized diffusion region is thought to originate from 
microscopic plasma processes occurring on scales many orders of magnitude smaller than 
the global scale. Thus, it is necessary to address the huge scale gap between the micro and 
the global scales while maintaining efficient global energy conversion. 


BiskampI fjl986[) argued that Sweet-Parker sheets with aspect rat io fi.e., length to thic 


k- 


ness) 


exceeding 100 are vulnerable to secondary tearing instability. 


Shibata and Tannma 


(1200 ll) developed a fractal reconnection model based on this concept. They reasoned that 


once plasmoids are formed by the primary tearing instability, plasmoid ejection and growth 
stretch the intervening current sheets, increasing the aspect ratio of the current sheets. 
Eventually these current sheets become unstable to secondary tearing, and disintegrate 
into chains of smaller-scale plasmoids and current sheets. By this stepwise process, an 
initially long and thick current sheet can reduce to the scale of the ion Larmor radius or 
ion inertial length. This hierarchical structure can plausibly couple the largest and smallest 
scales. High-resolu tion numerical simulations have proved the feasibility of this scheme 


(iBarta et ah 


20 nil) and also showed that the reconnection rate b ecomes independent of the 


Lnndquist number when S' > S'c ~ 10^ (e.g.. 


Loureiro et ah 


20121 ) . 
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This fractal reconnection model assnmes translational invariance in the direction 
perpendicnlar to the reconnection plane. This assnmption is violated in the presence 
of three-dimens i onal (3D) instabilities, snch as the kink-like instability reported by 


Dahlburg et ah 


(119921) ■ These anthors fonnd fast growth and saturation of the tearing mode 


along the anti-parallel magnetic held. Meanwhile, the oblique modes (with a component 
perpendicular to the tearing plane) continue growing and dominate in the later phase. 
Moreover, in the absence of a guide held, this kink-like instability will suppress the 
coalescence instability in the tearing layer, completely breaking the initial laminar structure 


flDahlburg and Einaudill2002[) . This instability is expected to ahect the fractal reconnection 


in 3D scenario. 


Another concern is the emergence of multiple tearing layers (IGaleev and Zelenvi 


1977h . 


When a sheared current sheet is subjected to a multi-modal perturbation (comprising modes 
ki, k2, ..., kn), multiple layers whose local magnetic held orientation satisfy k • B = 0 can 
emerge and are expected to interact when they become sufficiently close. Interaction should 


disrupt the 


Onofri et ah 


aminar structures of the layers, altering their internal fractal reconnections. 


(120041) studied the nonlinear evolution of several tearing layers coexisting 


inside a sheared current sheet in an incompressible plasma. They found that although 
the two-dimensional (2D) mode grows more rapidly than the emerging modes in a linear 
analysis, the emerging modes are oblique modes at an early stage. The resulting energy 
cascade gradually extends outwards from the current sheet center, leading to a hnal 


turbulen 

islands. 


state. Mor e over, the inverse energy transfer implies coalescence of the magnetic 


Landi et al. 


(120081) further conhrmed the three-dimensionality of the initially 


emerging modes in a compressible plasma simulation. Thus, the preferred source of 
plasmoid instability growth is change d and the reconnecti on e nhancement concep t should be 


modihed in the 3D case. But neither 


Onofri et al 


(2004) nor 


Landi et al. 


(120081) explicitly 


discussed the changes in the reconnection rate. For this purpose, we seek more details of 
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the magnetic energy consumption. 

This paper presents a detailed study of magnetic reconnection under tearing instability 
in a 3D resistive compressive MHD environment with a hnite guide field. In Section 2, we 
introduce the simulation model. In Section 3 and 4, simulation result and its discussion are 
presented. Section 5 makes a summary of the whole study. 

2. Simulation Model 

The 3D resistive MHD equations are solved in Cartesian coordinates. Viscosity, gravity 
and heat conduction are neglected for simplicity. The basic equations are as follows: 


^ + (v ■ V)p = -p(V ■ v) 

(1) 

i \ ^ J X B 

Vp+ ^ 

(2) 

J = X B 

An 

(3) 

^ + (v ■ V)p = -7P(V ■ v) + (7 - l)pj2 

(4) 

(9B 

—— = V X (v X B — cpJ) 
dt ^ ^ 

(5) 

P = -ksT. 

(6) 


Here, ks is the Boltzmann constant, rj is the resistivity, 7 = 5/3, c is the speed of light, and 
rh is the mean particle mass. For a fully ionized hydrogen gas, rh = 0.5mp, where rUp is the 
proton mass. The other variables take their usual meanings. 

All quantities are normalized by their characteristic values. The length scale is the 
width of the initial current sheet (= 5). The time scale is normalized by tA = 5 /vaqi 
where vao is the Alfven velocity. The initial mass density po normalizes the mass density. 
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The magnetic field is normalized by Bq = vao-s/Po- Consequently, the current density is 
normalized by Jo = cBq/S and the plasma pressure is normalized by po = -Sq- 



x/6 


Fig. 1.— Initial magnetic field plot along y = 0,z = 0. Solid line is the normalized 
By{x, y = 0, z = 0) and dashed line is the normalized Bz{x, y = 0, z = 0). 


The initial magnetic field comprises an anti-parallel component By{x) and a uniform 
finite guide field 5^, as plotted in Fig. [TJ 


B 





0y OtByQG^, 


(7) 


where Byo = \/4^i?o and a controls the magnitude of the guide field B^. Here we adopt 
a = 0.1 in the basic model. A magnetic field with non-zero a shears the structure along the 
a:-axis. The magnetic diffusivity is assumed to be spatially and temporally constant with a 
magnitude of = c^?7/(47r) ~ 3 x Resulting Lundquist number is thus defined 

by the Alfven speed va = uao\/1 + S = VALy/{2fj) ~ 2.6 x 10^. 


The simulation begins with a uniform mass density po throughout the simulation space. 
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Pressure balance between the gas and magnetic field is also maintained in the whole box, 
giving: 

P+|^ = ^(l + a")(l + /3) (8) 

in which [3 is the ratio between the plasma and magnetic field pressures at the minimum 
plasma pressure; we set /9 = 0.2. 


To initiate the reconnection, we add a random velocity perturbation with small 
amplitude {vx,Vy,Vz ^ 1 x over the whole simulation domain. Short-wavelength 

perturbations are eliminated. The remaining waves satisfy l/ky + l/kl ^ 1165^, also limited 
by I kyd 1^ 2 and | kzS 2. 


The simulation box size {L^ x Ly x Lz) is 105 x 245 x 65, containing 240 x 768 x 192 
grids. To resolve the tearing layer, we construct non-uniform grids with Ax ^ 0.025 along 
the x-direction. Uniform grids with Ay = Az = 0.031255 are constructed in the y- and 
^-directions. Periodic boundary con ditions are imposed on each border. The calculations 


are performed in CIP-MOCCT code flKudoh et ah 


(ILapidus 


19991 ) with an artificial Lapidus viscosity 


1967h . 


3. Simulation results 

The efficiency of the magnetic energy conversion is approximately five times higher 
in our 3D reconnection than in the Sweet-Parker reconnection. We propose a detailed 
mechanism of reconnection enhancement, referring to our new model as the “shock-evoking 
positive-feedback” model. 
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Fig. 2.— Current density |J|/Jo in the z = 0 plane (upper panel) and in the x = —0.35 
plane (lower panel). 


3.1. Global structure 


The upper and lower panels of Fig. [2] show the temporal evolution of the current 
density |J|/Jo in the a;|/-plane (z = 0) and the zy-plane (x = —0.35), respectively. In the 
early phase (Fig. [2](a)), diffusion dominates and no clear patte rn is observed ins ide the 
primary (central) current sheet. Gradually, tearing instability flFurth et al.l 119631) emerges 
and two chains of smaller diffusion regions (or reconnection sites) develop and grow adjacent 
to a; = ±0.35 (Fig. [2](b)). Within these chains, the current density is much higher than the 
surroundings. The chains form a web-like pattern across the global current sheet. Fig. [3] 
and Fig. m present the plasma density, pressure, velocity, magnetic held, and current density 
at t = 400tyi and 700tA respectively, in the z = 0 plane. The A notation dehnes a difference; 
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Fig. 3.— Global patterns of variables on the z = 0 plane at t = 400t^; plasma density 
difference (Ap), pressnre difference (Ap), velocity components {vx,Vy,Vz), magnetic field 
components and differences {ABy, AB^), current density component {Jy, J^) 

for example, Ap is defined as Ap = p — PiniUai, where Piniuai is the initial value of p. The 
structures are well-organized in each plot of Fig. [3l but becomes distorted by flows in the 
later phase (Fig. 0]). 

Oblique lines (diffusion lines) are formed by the diffusion regions in the zp-plane, as 
shown in the lower panel of Fig. |2](b). The positions of these diffusion lines along the 
x-axis are determined by resonance layers (resonant tearing layers), which arise from the 
periodic boundary conditions in the y- and 2 :- directions. This slab structure resembles 
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Fig. 4.— Global patterns on the z = 0 plane at t = 700t^. Variables are defined in the same 
way as Fig. [31 


the cylindrically symmetric rational snrface in a tokamak (e.g., 


Behan 


2006 1. The safety 


factor q, which denotes the period of the magnetic field lines cycling (sheet-wise) across the 
z|/-plane, is calcnlated along the x-axis as follows: 


q{x) 


Ly {Bz)x 
{By') X 


(9) 


Here, Ly and are the box sizes in the y- and ^-directions, respectively. {By)x and {Bz)x 
are the mean magnetic field components on the zj/-plane at the corresponding point x. 






































11 


Strictly speaking, provided that the q value in a zy-]Aajie satishes 

m 

q = — 

n 

where m and n are integers, the wave signal is strengthened and the layer develops into 
a resonance layer. Clearly numerous layers are expected. However, due to the different 
growth rates of these modes, a limited number of layers can resonate locally. When m and n 
are large, the corresponding wavelength is small; hence, the instability is easily suppressed 
by magnetic tension force. 



The safety factors calculated by Eq. (9) at t = 250^^ are plotted as black solid line 
in the upper panel of Fig. |5l To find the most unstable mode in the box, the Fourier 
transformation of the current density |J|/Jo on each z|/-plane is calculated along the a:- axis. 
The signal intensity is calculated as 



Here, ky and kz are wavenumbers (defined as A; = 1/A, where A is the wavelength) in the 
y- and z-directions, respectively. The intensity increases and fades when approaching and 
departing a resonance layer, respectively. The changes in /j and |J|/Jo along the x-axis are 
demonstrated in five slices at a; = —0.55, —0.35, 0, 0.35, 0.55 (see middle and lower panels of 
Fig. [5]). The black dash-dotted lines in the middle panels indicate the vectors perpendicular 
to the local magnetic field. On certain resonance layer [x = ±0.35), the line coincides with 
the maximum signal since it must satisfy k ■ B = 0, where k = kyey ± kz^x- As By changes 
sign across x = 0, the dominant k in the resonance layers reverses the sign of one of its 
component while preserving another. The blue solid line in the upper panel of Fig. [5] plots 
the maximum of /j on each z^z-plane along x-direction, max^^ ^^^(/j), as a function of x. We 
hnd that maxk^^kyifi) peaks several times inside the primary current sheet. Comparisons 
with the q curve reveal that these peaks correspond to resonance layers with q = 1 and 
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corresponding layers. 

g = 2 (red and orange dashed lines respectively in npper panel of Fig. |5]). Resonance layers 
with other q valnes are missing in this model, possibly becanse that the large-g resonance 
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layers near the center {x = 0 ) are separated by less than the spatial resolntion of the 
simnlation. It is also likely that the most nnstable tearing mode in this conhguration is the 
3D mode {kz 7 ^ 0) rather than a 2D mode [kz = 0), and the modes o n layers with ^ 3 


grow at smaller rates than the modes on layers with g = 1 and g = 2 flBaalrnd et alJ 


201211 ■ 



selected magnetic 
field lines 


\J\/k 


0 

y/S 


- 36 


Fig. 6 .— Coherent structnre of flnx tubes extracted from the isosurfaces of current density 
|J|/Jo = 0.12 at f = SOOfyi. The plotting area in the y- and ^-directions is expanded by a 
factor of 3. Translucent white surface represents the a; = 0 plane. Colored solid lines are 
selected magnetic held lines. 

As the guide held is uniform, no null points or sheets exist. The magnetic held lines 
reconnect across the tearing layers and only part of the held strength is consumed. These 
component reconnections lead to twisting reconnected held lines that cross the resonance 
layer to and fro in 3D space. The held lines form coherent hux tubes whose width along 
the z|/-plane approximately equals the wavelength A of the tearing mode (A = 2.96 and 
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5.45 on q = 1 and q = 2 planes, respectively). The flnx tnbe strnctnre can be extracted 
from the isosnrfaces of the cnrrent density, as shown in Fig. [H] (since periodic bonndary 
conditions are assnmed, the plotting area in the y- and 2 ;-directions is extended by a factor 
of 3). Twisting held lines wind aronnd the hnx tnbes. The tilting angle of the hnx tnbes 
reverses across the a; = 0 plane (indicated by the translucent white surface in Fig. [6]). 

Later, as the local reconnection events gradually develop in the individual diffusion 
regions, the hnx tubes thicken along the x-direction. The outhows from the diffusion regions 
are sufficiently accelerated to distort the well-organized structure (Fig. [2](c)). The hnx 
tubes centered at the negative-a; side collide and coalesce (Fig. 121(d)). Finally, only spatially 
extended current sheets and large hnx tubes remain (Fig. [2](e)). 

Actually there are three current sheets coexisting inside the whole simulation box. It 
could be seen clearly that the current sheets at the boundary have no signihcant influence 
on the current sheet at the center (Fig. [7j). Furthermore, We test a reconnection model with 
twice the length along a:-direction to see the effect of the boundary and the result does not 
change. 



Fig. 7.— Current density |J|/Jo and velocity Vx/vao in the z = 0 plane by comparing plot 
ranges from —35 to 35 and plot ranges from —55 to 55. 









15 


3.2. Reconnection rate 


To understand the reconnection efficiency, we acquire the reconnection rate calculated 
as follows: 


Ma = 


d 

dt 


err^dV 


2LyL^-^VAo 

dvr 


( 12 ) 


In the numerator of Eq. (12), the magnetic energy density em is integrated over the volume 
inside the primary global current sheet to obtain the total magnetic energy. The boundary 
of the primary current sheet is determined from the specihed critical plasma pressure 
Pc ~ O.OOSRg. The energy inside the considered sheet is conserved by adding the surface 
flux term along the a:-direction in Eq. (4). The result is plotted in Fig. [HI The reconnection 
rate is rapidly enhanced after t = 450^^- 



Fig. 8.— Reconnection rate in the basic model (black solid line), whose maximum is about 
5 times of the Sweet-Parker type reconnection simulated in the same 3D box (not plotted). 
2.5D simulation with the guide held (green dash-dotted line), 3D simulation without a guide 
held (orange dash-dotted line), and 3D simulation with a single mode (red dash-dotted line). 
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To check whether 3D reconnection is sufficiently fast, we simulated the Sweet-Parker 
type reconnection in the same 3D box. This model begins with four large vortices 
imitating the inflow-outflow coupling in a single reconnection. The velocity components are 
constrained by v^iVy ^ 2 x and Vz = 0. Translational invariance is maintained 

at the beginning. The Sweet-Parker reconnection rate, which is not plotted in Fig. [HI is 
0.00026, which is approximately 1/5 the reconnection rate of the basic model. 

To better quantify the effect of the third dimension on the reconnection, we conduct 
a 2.5D simulation with a guide held (green dash-dotted line in Fig. [H]). The initial velocity 
held is the initial velocity held on the z = 0 plane in the basic model. Because kz = 0 (or 
Lz oo) throughout the simulation box, the resonance layer is limited in a: = 0 plane. A 
secondary instability (plasmoid instability) starts to develop at the end of this simulation, 
and the reconnection rate increases later than in the basic model. In comparison, no hner 
hlamentary structures are detected in the basic model. Then the reconnection rate increase 
in the basic model is due to a secondary instability other than plasmoid instability. 

To generate oblique tearing layers, we also require a hnite guide held. If the initial 
magnetic held has no z-component {a = 0), the tearing mode can grow only along the 
a: = 0 plane, across which the magnetic held reverses. The result of a 3D simulation without 
the guide held is plotted as the orange dash-dotted line in Fig. [HI In this model, the initial 
velocity is randomly perturbed as in the basic model. Clearly, the reconnection is much 
lower in this simulation than in the basic model. 

As shown in the previous result, the basic model admits several tearing layers coexisting 
in the same current sheet. To understand the importance of these multiple layers, we 
conduct a 3D simulation with a single mode (k) in the perturbed velocity held (red 
dash-dotted line in Fig. [H]). The selected k is the most unstable mode in the leftmost 
resonance layer (at the negative x-side) of the basic model. The reconnection rate in this 
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Fig. 9.— Development of total consumed magnetic energy (blue solid line) and total kinetic 
energy (green solid line). Red dash-dotted line represents exponential growth. 

model is ~ 20% that of the basic model throughout the same development period. The 
reconnection within each diffusion region gradually saturates rather than largely increases. 
This result conhrms that the reconnection rate is enhanced by multiple resonance layers 
coexisting in one current sheet. 

To quantitatively examine the energy conversion efficiency, we plot the temporal 
evolutions of the normalized total kinetic energy Ek/B^, and the consumed total magnetic 
energy \AEm\/B q of the central global current sheet in Fig. [Hi The kinetic energy grows 
exponentially from t = 200^^ to t = 450t^. From exponential htting, we hnd that 
Ig^A ~ 0.039, where 'jg = d{lYLEi.)/dt. Comparing these results with the reconnection rate, 
it is easily seen that rapid enhancement immediately follows the exponential growth of 
kinetic energy (at t = 450t/i). Moreover, the tearing layers are already recognizable around 
this time (see Fig. [2](b)). Together with the above hndings, this implies that the interaction 
between fully grown tearing layers triggers the fast consumption of magnetic energy. 
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Regarding the different reconnection rates plotted in Fig. [HI we conclude that the 
interaction between different resonance layers is crucial for the enhanced reconnection in 
later phase. 


3.3. Positive-feedback system 


On closer examination, the diffusion regions on multiple tearing layers are observed 
to form a web-like pattern across the current sheet. Fig. [TU] plots the current density 
|J|/Jo on the z = 0 plane at various times of the simulation {t = 450t^, t = 500t^, and 
t = 550t^), overlaid with the velocity flows (white arrows) and the ^-component of the 
vorticity oJz/ojq, where ojq = vao/^- On this plane, the small diffusion regions form an 
asymmetric structure with a zigzag pattern. Within this structure, the diffusion regions 
on different resonance layers are apart with each other in ^/-direction. Examining the local 
stream, we hnd that the outflow from one reconnection site diverts and feeds into the inflow 
region of the reconnection site on a different resonance layer. Meanwhile, a portion of the 
reconnected magnetic flux is transported, where it can again participate in reconnection. 
The transportation is regulated by the outflow from the diffusion region. The outflow 
strengthens as the local reconnection proceeds, implying that faster flux transportation 
correspondingly enhances the inflow, thus accelerating the local reconnection. This coupling, 
named as positive-feedback system, is generated by the existence of multiple tearing layers. 


in contrast to the se 
reported by 


f -feedin g system inside a current sheet with a single reconnection layer 


Lapental (120081) . The coupling of inflow and outflow regions across the current 


sheet can be identified by the coalescence of branch-like structures, which extend from 
the individual diffusion regions in the contour plots of An example of such coupling 
is delineated by the black dashed rectangle in the lower panels of Fig. [TOl Although this 
feature is gradually deformed by the turbulence developing inside the current sheet, the 
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Fig. 10. — Snapshots of current density |J|/Jo and the ; 2 -coniponent of vorticity oJz/ojq on 
z = 0 plane at / = 450/^, 500/^, and 550/^. White arrows represent the flow patterns. The 
vector scale is shown at the bottom of the left upper panel. 

secondary-transportation of magnetic flux is maintained, ensuring that reconnection can 
proceed. 

As shown in Fig. |5l if of the most unstable modes retain their sign, ky of these 
modes change sign at opposite side of the current sheet center. This implies that the tilting 
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Fig. 11.— Upper panel: 3D surface of current density |J|/Jo = 0.12 in the simulation box. 
Black translucent surfaces are the planes at z = —1.3755, —0.68755, and 0 at t = 500^^- 
Middle panel: Snapshots of |J|/Jo on the planes shown in the upper panel. White arrows 
represent the flow patterns. Lower panel: Corresponding Uz/oJo contour plots. 
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Fig. 12.— Cartoon of the layout of diffusion regions across the global current sheet. Red 
and purple arrows represent the characteristic flows on the xy-plane. 


angles of the diffusion lines (dehned as the angles between the orientations of the diffusion 


lines and the xz-plane) differ at either side of the global sheet, as shown in the flux tubes 
of Fig. [6l The tilting angle is dehned as: 



(13) 


Therefore, the conhguration of the diffusion regions changes in different xy-planes along 
2;-axis. Fig. [TT] plots the same variables as Fig. [10] but at different ^(-positions {t = 500fA)- 
Unlike their appearance in the z = 0 plane, the diffusion regions on the different tearing 
layers shift in the y-direction and sometimes show a symmetric structure (in which diffusion 
regions on different resonance layers align along the a:-direction). In this conhguration, 
the two dihusion regions do not cooperatively interact; thus no feedback character is 
established. The transition can be tracked by observing the pattern delineated by the 
dashed rectangle in the lower panel of Fig. [TTl To an observer moving along the z-direction. 
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the diffusion regions form an alternating asymmetric-symmetric-asymmetric structure. A 
cartoon of these structures are presented in Fig. O 

To examine the global structure of the positive-feedback system built by multiple 
tearing layers, we plot the 3D and 2D spatial distributions of the local Vx (see Fig. [T3H . 
Once the outflow region of a certain reconnection site couples with the inflow region of 
another site, the surface of relatively high Vx continuously extends across the current sheet 
center. Here, we select a cutoff surface of Vx/vao = ±0.0005, approximately 20% of the 
maximum absolute value of Vx/vao throughout the box. The dashed rectangles in panels 
(a) and (b) of Fig. [12] highlight where the surface crosses the black translucent plane (the 
central plane of the current sheet at a: = 0). On this plane, the Vx contour plot exhibits an 
oblique checkered pattern, as it is presented in Fig. [T3fdb The same features are highlighted 
in the same regions. Noted that Vx is not enhanced at the corners of the checkered pattern, 
because the diffusion regions are symmetrized at those points. Where the diffusion regions 
become asymmetric, the positive-feedback system establishes provided that a little shift 
occurs between the diffusion regions. Thus, the Vx is distinctly enhanced along the borders 
of the checks. 

In summary, the diffusion regions in different resonance layers couple via their inflow 
and outflow regions. This conhguration, the positive-feedback system, globally occurs inside 
the box, and underlies the rapid reconnection observed in our model. 


3.4. Global inflow and slow-mode shocks 

The enhanced reconnection rate in the basic model implies that a large amount of 
magnetic energy is transported into the current sheet. To evaluate the rate of energy 
transport, we investigate the global inflow along the x-boundary of the central current 
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I -0.0005 
■ 0.0005 


Fig. 13.— (a)-(c) Surface plots of and | J| in the 3D box. Grey surface represents | J|/Jo = 
0.12. Blue and red surfaces denote v^/vao = —0.0005 and 0.0005 respectively. Black, 
green and yellow dashed rectangles highlight the coupling of inflow and outflow regions, (c) 
combines panels (a) and (b). (d) is the 2D contour plot of Vx on the translucent plane in 
(a)-(c). The same features are highlighted in the the same regions. 
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<|J|/Jo), onz = 0 



Fig. 14.— (a) Normalized Poynting flux at 700tA on the z = 0 plane. Black contours 
indicate the enhanced current density |J|/Jo ^ 0.8. Current sheet boundary (pc = O.OOSBq) 
is marked by the red and blue dashed lines to the left and right of a: = 0, respectively, (b) 
Averaged global inflow (black solid line) and reconnection rate (black dash-dotted line), (c) 
and (d) Spatially averaged normalized (E x B)^; and |J| inside the current sheet indicated 
by the dotted double arrows in (a). 

sheet (Fig. [T^ . Fig. [TU(a) exhibits snapshot of the normalized Poynting pattern (E x B)^ 
at the time of Fig. [2](d) (namely, at f = 700tyi)- The black contours indicate where the 
current density is locally enhanced. The surfaces assumed in the average inflow calculation 
are highlighted by the colored dashed lines to the left and right of the global current sheet. 
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They are the same surfaces used in the reconnection rate calculation. The average global 
inflow is determined as follows: 

^(E X B) -dS 


Vgx — 


‘^B^LyLz 


(14) 


where E is the electric held and dS denotes the surface on global current sheet boundary. 
We consider the dot product to be only the x-component of E x B, thus |dS| = dydz. 
This quantity rehects the transport rate of the Poynting hux into the current sheet in the 
x-direction. As shown in Fig. [TTf b). the Vgx changes smoothly over time. Panels (c) and (d) 
of Fig, im present time prohles of the normalized E x B and |J|, respectively. Both variables 
are spatially averaged inside the current sheet as follows (where / denotes a variable) 


(/).= 


(15) 


Xr{y)-xi{y) 

In Eq. (15), Xi{y) and Xr{y) denote the positions of the left and right boundaries, 
respectively, of the current sheet at (?/, z = |/,0). The diffusion region is asymmetric in 
the z = 0 plane at f = 700^^, as observed in Fig. ITTfa). Consequently, the inflow regions 
extending from the current sheet are also asymmetric. This pattern becomes increasingly 
obvious as local reconnection is accelerated by positive-feedback system iFig. [m c)b When 
entering the turbulent state (after ~ f = 600^^; see Fig. imd)). the dense regions essentially 
overlap with the structures of Fig. im a) and (c), and the current density inside the current 
sheet is universally enhanced. Although the inflow structure is segmented, the global effect 
can be regarded as a single long diffusion region. Moreover, beyond t = 400^^, the Vgx 
well correlates with the reconnection rate (black dash-dotted line in Fig. fMf b)). The 
enhanced Vgx suggests that, on macroscopic scale, magnetic energy is carried into the global 
current sheet by the flow and converted into other forms of energy. This suggestion is 
supported by the rising reconnection rate in later phase. The average global inflow exceeds 
the reconnection rate at approximately t = 700f^, indicating that a small portion of the 
magnetic energy is stored before being gradually consumed. 
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x/5 x/5 


Fig. 15.— Plots of normalized current density, plasma pressure, magnetic pressure, entropy, 
plasma density and velocity across {y^z = 0,0), showing the development of a slow-mode 
shock. Variables are plotted at t = 640^^ (black solid lines), 670tyi (blue solid lines), TOOt^ 
(orange solid lines). Crosses indicate the grid points. Red dash-dotted line is the detected 
shock front at t = 700tA- Black horizontal dotted lines denote the shock speed in lab-frame. 


When the inflow enhancement is high, slow-mode shocks gradually form between the 
inward flow and the outwardly growing flux tube. Along with reconnection in the diffusion 
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region, these shocks are a likely mechanism of magnetic energy conversion. Shocks are 
identihed by comparing the npstream and downstream qnantities in the shock frames at 
selected points. Slow-mode shocks in the simnlation domain are reqnired to satisfy the 
following criteria: 


1. Rankine-Hngoniot relations: The Rankine-Hugoniot relations in 


not deviate by more than 30% of their npstream valnes (jSaito et ah 


he downstream must 


19950; 


2. Velocities must satisfy vau ^ Vnu ^ Vsiu, where Vnu is the normal component 
(perpendicular to the shock front) of the upstream velocity, vau and Vsiu are the Alfven 
speed and the slow-mode wave phase speed, respectively, in the upstream; 

3. Vnd ^ Vsid, where Vnd denotes the normal component of the downstream flow speed and 
Vsid is the slow-mode wave phase speed in the downstream; 

4. The B and v of the upstream and downstream must be co-planar (within 10°). 


Fig. [15] shows the development of a slow-mode shock wave. Plotted are the current 
density Jz/Jo, plasma pressure p/po, magnetic pressure pm/po, entropy (s — so)/cy, plasma 
density p/po, plasma velocity Vx/vao, Vy/vAo and Vz/vao across (p, z = 0,0) in laboratory 
frame at three time points (for tracking the transition). The entropy is normalized by 
So = cyln(po/Po), where cy is the specihc heat at constant volume. The shock speed in 
the laboratory frame is indicated by the black dotted line in the plasma velocity plots. 

At t = 640f/i (black solid lines in Fig. [T5|), the density is compressed around x = —6. 
The local Jz is minimized at the x-axial boundary of the flux tube. At this moment, the 
local x-directional flow on either side of the flux tube boundary is relatively weak. As the 
reconnection becomes more efficient, the flux tube thickens and expands in the negative 
x-direction. The relative flow between inside and outside of the tube boundary increases, 
further compressing the local plasma (compression occurs around x = —1.35). Eventually 
(at t = 700f/i), the shock criteria are reached. The approximate position of the shock front 
is indicated as red dash-dotted line in Fig. [151 We fail to hnd the rotational discontinuity 
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previously reported in reconnection with a guide field (e.g., 


Longcope et al. 


200911 ■ 


Normalized Poynting 
flux reduction 



t/tA 

(a) 


slow-mode shock 
coverage percentage 



t/tA 

(b) 


Fig. 16.— Temporal evolutions of (a) normalized energy conversion capability of shock Ms 
between the upstream and downstream of slow-mode shocks throughout the simulation box, 
and (b) percentage area coverage of slow-mode shock. 



Fig. 17.— All detected shocks (indicated by red dots) are located approximately by down¬ 
stream position at t = 700f^. Blue surface represents = 0. 
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Along with local Ohmic heating, these slow-mode shocks are considered as an 
additional source of magnetic energy consumption, as shown in Petschek’s work. The 
energy conversion capability of shock (reduction of Poynting flux across the shock) is 
calculated as: 

E[(ExB)|„-(ExB)U.dS 

=- wmr. -■ 

The hrst and second terms of numerator denote the upstream and downstream Poynting 
vectors, respectively. The Mg is plotted as a function of time in Fig. fTbf a). The efficiency 
of these shocks is greatly increased at the end of the simulation. To understand the global 
layout of the slow-mode shocks, we calculate the total percentage area coverage SC of the 
shock: 


SC = X 100% (17) 

4:LyLz 

where Sgs is the total area covered by the shock. The coefficient “4” in the denominator 
indicates that four tearing layers coexist in the central current sheet. The temporal 
evolution of SC is plotted in Fig. [TOfbh Both Ms and SC increase from t = 450fyi. 
Especially in the later phase, when the global inflow booms (after t = OSOf^), large extent of 
slow-mode shocks with efficient energy conversion capability are detected. The approximate 
downstream positions of all shocks found in the domain are highlighted in red in Fig. [T71 
The blue isosurface is the flux tube surface with = 0. 


In summary, multiple tearing layers construct a positive-feedback system that promotes 
reconnection. Large amount of magnetic energy is transported into the current sheet and 
converted to other forms of energy. Outwardly growing flux tubes collide with the strong 
inflow, gradually arousing slow-mode shocks along the tube boundaries. These shocks 
further enhance the reconnection. Therefore, we refer to our model as the “shock-evoking 
positive-feedback” model. 
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4. Discussion 


In our present model, the reconnection rate is enhanced by non-linear interactions 
among multiple (tearing) layers. Our observed int eractions resemble the double tearing 


mode (DTM) reported in the previous studies (e.g., 


PurthetaL 


19731) . Several comparisons 


are worth a brief mention here. First, the unperturbed state of DTM has multiple 
independent current sheets, whereas our model begins with a single current sheet. The 
multiple current layers are consequent to the growth of the tearing layers by mixed 
perturbations. Such growth of multiple layers requires a 3D system and a moderate guide 
held. Second, DTM is essentially a 2D process because translational invariance is assumed 


along the guide held direction (e.g.. 


Janvier et ah 


201 ll) . Therefore, the arrangement of 


dihusion regions are maintained in that direction. In our model, the arrangement of dihusion 
regions dihers among z-positions (see Section 3.3) and periodically changes from symmetric 
to asymmetric. Because of this non-uniformity, positive-feedback system is intermittently 
distributed. Our model is expected to consume magnetic energy less efficiently than DTM. 
Nevertheless, our model substantially enhances the reconnection rate. Third, as also shown 
in Section 3.3, the energy release may be slow in regions of symmetric arrangement of 
dihusion regions from diherent tearing layers. Notably, howeve r, reconnection proceeds in 


symmetric DTM structures albeit very slowly (lYan et ah 
with our model. 


1994) . This result is consistent 


In applying our model to realistic systems, such as solar hares, we must reconsider 
the boundary ehects and the resistivity, which is much larger in our model than in solar 
phenomena. In our simulation, the selection of the tearing layers is highly inhuenced by 
the box size; in other words, by the aspect ratio of the initial current sheet. In a much 
larger system enclosing a current sheet with an extremely large aspect ratio, the evolution 
should be inhuenced by the boundary ehects. If the plasma (3 is high in the current sheet. 
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local Alfven speed is slow. The evolution may be locally determined and the system will 
behave as an open system with free outflows. In low-/! plasma, the evolution is globally 
determined and influenced by the boundary conditions at both ends of the held lines. 
Meanwhile, the most unstable tearing mode depends on both the resistivity and the guide 
held. However, we argue that once the emerging multiple tearing layers are sufficiently 
close to interact (separated within several 1/|A'|, where 1/|A'| is the width of tearing 
layer), a positive-feedback system is built. Conversely, if the multiple layers approach 
very closely (within 1/|A'|), they will behave as a single layer. These propositions require 
investigation in further study; for instance, we should elucidate the scaling law that relates 
the Lundquist number to the reconnection rate, survey the guide held strength parameter 
and investigate the boundary condition. We must also improve the resolution and precision 
of our simulation. 


Finally, we propose an extended model based on our “shock-evoking positive-feedback” 
model, which operates similar to 2D hierarchical reconnection. In the later phase of our 
simulation, long and thin current sheets with sheared ma gnetic structures remai n near the 


center, which are candidate structures for further tearing. 


Daughton et ah 


f 201lh simulated 


reconnection in a single current sheet under a guide held ehect in the kinetic regime. They 
reported secondary hux ropes resulting from tearing instability. No secondary hlamentary 
structure is observed in our model, likely because of the low resolution. Higher-ranking 
hlaments are expected in our model when the Lundquist number and resolution are 
improved. Once the extended current sheet becomes unstable to the tearing instability, 
there are two possibilities. If the most unstable mode is centralized, a single tearing 


layer can be identihed, a s occurs in 2D 


unstable mode is oblique flBaalrud et ah 


plasm oid instability. Conversely, if the most 


2012n . multiple tearing layers should emerge. The 


coexisting multiple layers resemble the primary current sheet structure. Once established, 
the positive-feedback system enhances local reconnection. Flux tubes at the same side 
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Fig. 18.— Illustration of hierarchical structure in the “shock-evoking positive-feedback” 
model. 


grow and merge with each other. Slow-mode shocks are evoked along the outer boundary, 
while extended current sheets form between the coalesced flux tubes. These remaining 
current sheets might be subject to further tearing instability. Therefore, in a series of steps, 
the global structural scale of the diffusion region could reduce to microscopic size (such 
as the ion inertial scale) by alternate applications of 2D plasmoid instability and the 3D 
“shock-evoking positive-feedback” model. We regard this fractal structure of the current 
sheet generated by the “shock-evoking positive-feedback” model as an extension of the 
hierarchical structure of 2D plasmoid instability. Our scenario is illustrated in Fig. [THl 
assuming that multiple tearing layers are realized in each rank. 












5. Summary 


In this study, we simulate a 3D MHD magnetic reconnection across a symmetric 
current sheet with a hnite guide held. We assume a uniform resistive environment and 
randomly perturb the initial velocity held. When relaxing the variation along the direction 
of the guide held, we hud that the guide held increases the reconnection rate of the whole 
current sheet by several times, relative to the simulation with no guide held. Most unstable 
tearing modes, with components in the sheet-wise direction, emerge on multiple resonance 
layers. The dihusion regions on these layers establish a zigzag pattern that couple the 
inhow region and outhow region on diherent layers. Consequently, the outhow is diverted 
into the inhow region, ensuring that reconnection proceeds in the opposite layer. Gradually, 
the inhow from outside of the global current sheet also becomes accelerated by continuous 
activation of individual reconnection site. This enhanced inhow arouses slow-mode shocks 
along the outer boundary of the current sheet, further promoting energy conversion. We 
refer to our model as the “shock-evoking positive-feedback” model. 

The Lundquist number is much smaller in our numerical experiment than in the 
solar corona. Therefore, our model is not directly applicable to real solar activities. 

To understand the feasibility of the “shock-evoking positive-feedback” in low resistive 
environments, we must conduct a parameter survey on the dihusivity magnitude. Altering 
the resistivity and guide held strength would induce diherent most unstable tearing modes, 
thus changing the topology of the positive-feedback system. In future study, we will further 
evaluate the ehectiveness of our model by varying and f/. 

The limited resolution of our model precludes the detection of hner hlamentary 
structure. If the remaining long, thin current sheet becomes vulnerable to tearing 
instability, the internal sheared structure might develop into a higher-ranking positive- 
feedback system. The mechanism of the “shock-evoking positive-feedback” model would 
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then enhance the local reconnection rate. In this way, the pattern could be repeated on 
ever smaller scales, eventually reaching the microscopic scale in environments of large 
Lundquist number. Therefore, we have potentially expanded 2D fractal reconnection into a 
3D hierarchical structure with a high energy conversion rate. This capability of our model 
needs to be tested on very hne simulation grids, another goal of our future work. 

This research was conducted using the Fujitsu PRIMEHPC FXIO System (Oakleaf- 
FX,Oakbridge) in the Information Technology Center, The University of Tokyo. Numerical 
computations were in part carried out on Cray XC30 at Center for Computational 
Astrophysics, National Astronomical Observatory of Japan. This research is supported by 
Leading Graduate Course for Frontiers of Mathematical Sciences and Physics (FMSP) and 
JSPS KAKENHI Grant Number 15H03640. 
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